%{
AUTHOR: Felipe Arteaga
-------------------------------------------------------------------------
PROJECT: Warnings
-------------------------------------------------------------------------
DESCRIPTION:
=========================================================================
%}


clearvars -except projectDir projectDirData fromMainWarningsPaper
clc;close all;fclose('all');feature('DefaultCharacterSet','UTF-8');

if(not(exist('projectDir','var')==1&&exist('projectDirData','var')==1&&exist('fromMainWarningsPaper','var')==1&&fromMainWarningsPaper))
    pcName=char(java.lang.System.getProperty('user.name'));
    if(strcmp(pcName,'felipe'))
        % PC Felipe
        myDir='/Users/felipe/Dropbox/';
        projectDir=[myDir,'git/warnings/'];
        projectDirData=[myDir,'projects/warnings/'];
        addpath(genpath([myDir,'/myMatlabFunctions/']));
    end
end
compileLatexTable=false;


%%
dirPlots=[projectDir,'/paper/figuresCL/playRDs/'];

% warning('Alternative dir plots')
% dirPlots=[projectDir,'/paper/figuresCL/RDs/Aux/'];

dirTable=[projectDir,'/paper/tablesCL/'];
dirData=[projectDirData,'/data/chile/'];

tableBandwidthComparison=false;
plotRDs=false;
titlePlotRD=true;
plotWithZoom=false;
savePlotRDs=true;
bandwidthTable='local'; % 'full','local','localAuto';

switch bandwidthTable
    case 'full'
        bwnote='Full bandwidth was used with 2nd order (local) polynomial fit. ';
    case 'local'
        bwnote='Bandwidth was set to +-0.1 with default local polynomial fit (1st order). ';
    case 'localAuto'
        bwnote='Bandwidth was set to automatic (bwselect with default options) with default polynomial fit (1st order). ';
end

%% Load data:
anho=2017; %0: pooled
if(anho==0)
    anhoStr='';
else
    anhoStr=sprintf('%i',anho);
end


notIn2020={}; % {'assignedToPref','acceptOffer','declineOffer','preferenciaAsign_1','placedInAnyPreferenceIniEnd','placedInAddedIniEnd','placedInAddedNoWaitlistIniEnd','placedInOldIniEnd','enrolledInAssigned','enrolledInAddedIniEnd','enrolledInAddedNoWaitlistIniEnd','enrolledInAddedNoRiskIniEnd'};

    load([dirData,anhoStr,'/inputRD'],'dataRD')


dataRD.rural=dataRD.newMarket<0;

% As asignment at .3 is control, just move (if any) in .3 to epsilon to the
% left:
dataRD.riskPopup(dataRD.riskPopup==.3)=.29999999;
dataRD.riskPopup(dataRD.riskPopup==.5)=.49999999;
dataRD.riskPopup(dataRD.riskPopup==.7)=.69999999;

% Avoid mass points at extremes:
dataRD.restrAll=dataRD.pobPopup==1&dataRD.riskPopup>.01&dataRD.riskPopup<.99;

% No need to be more specific, and put "Predicted placement risk of 1st
% attempt"
dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
dataRD.Properties.VariableDescriptions{'treatedPopup'}='Treated pop-up in first attempt';


if(false)
    %% % Popups of total population
    n_apps=[274990,483070,454415];
    
    fprintf('2018 In data: %.1f, in popup population: %.1f\n',height(d1.dataRD)/n_apps(1)*100,sum(d1.dataRD.pobPopup)/n_apps(1)*100);
    fprintf('2019 In data: %.1f, in popup population: %.1f\n',height(d2.dataRD)/n_apps(2)*100,sum(d2.dataRD.pobPopup)/n_apps(2)*100);
    fprintf('2020 In data: %.1f, in popup population: %.1f\n',height(d3.dataRD)/n_apps(3)*100,sum(d3.dataRD.pobPopup)/n_apps(3)*100);
    
    fprintf('\n ALL  In data: %.1f, in popup population: %.1f\n',height(dataRD)/sum(n_apps)*100,sum(dataRD.pobPopup)/sum(n_apps)*100);
    
end





%% RD

% Defines subpopulation in which I want to calculate RDs for outcomes
% defined soon.

assert(all(dataRD.cutoffAPI_ini(dataRD.restrAll)>0))

    

subpobs=cell(1,8);
s=1;
subpobs(s,:)={'Pooled',    dataRD.restrAll,'pooled','addAnyIniEnd','0.0','.28','.28','riskPopupCentered'};s=s+1;
subpobs(s,:)={'0.3',   dataRD.restrAll&abs(dataRD.cutoff-0.3)<.001,'c.3','','0.3','.28','.68','riskPopup'};s=s+1;
subpobs(s,:)={'0.5',      dataRD.restrAll&abs(dataRD.cutoff-0.5)<.001,'c.5','','0.5','.48','.48','riskPopup'};s=s+1;
subpobs(s,:)={'0.7',    dataRD.restrAll&abs(dataRD.cutoff-0.7)<.001,'c.7','','0.7','.68','.28','riskPopup'};s=s+1;
%posAll=1;




% Definition of outcomes for RDs. 3rd row defines the beginning of a panel
% in the table
% 4th row

if(ismember('DRiskLExpostIniEnd',dataRD.Properties.VariableNames))
    error('Borrar esto');
else
    dataRD.DRiskExpostIniEnd=dataRD.riskExpostEnd-dataRD.riskExpostIni;
end

% New: enrolled conditional on placed (short name because of stata)
dataRD.enrolledInAssignedCOA=double(dataRD.enrolledInAssigned);
dataRD.enrolledInAssignedCOA(dataRD.assignedToPref==0)=nan;
dataRD.placedInAnyPreferenceNWIE=dataRD.placedInAnyPreferenceNoWaitlistIniEnd;
dataRD.placedInAnyPreferenceNoRiskIE=dataRD.placedInAnyPreferenceNoRiskIniEnd;

%notIn2020=[notIn2020,'enrolledInAssignedCOA'];

% 5th column defines if IV is calculated or not.
outcomes={'changeAnyIniEnd','Any modification','cha','',0;...
    'addAnyIniEnd','Add any','aa','',0;...
    %'deltaProbPlacedNoRiskIniEnd','$\Delta$ prob. placed to zero risk','dpnr','',1;...
    %'deltaProbPlacedNoWaitlistsIniEnd','$\Delta$ prob. placed to uncongested','dpw','',1;...
         'schoolsAddedIniEnd','Schools Added','sa','',1;...
    %     'riskLastDayIni','True Risk initial attempt','ri','';...
    %     'riskLastDayEnd','True Risk final attempt','rf','';...
         'DRiskExpostIniEnd','$\Delta$ Risk','drEP','',1;...
    %     'riskExpostIni','True Risk initial attempt','riEP','';...
    %     'riskExpostEnd','True Risk final attempt','rfEF','';...
    %     'DRiskExpostIniEnd','$\Delta$ Risk','drEP','',1;...
    %     'numberOfNewUncongestedIniEnd','Uncongested schools Added','saUncong','';...
    %     'addAsFirstIniEnd','Add as first','af','',1;...
    %     'addInBetweenIniEnd','Add to middle','ab','',1;...
    %     'addAsLastIniEnd','Add as last','al','',1;...
    %     'deleteAnyIniEnd','Drop any','da','',1;...
    %     'changeOrigOrderIniEnd','Re-order','ro','',1;...
    %     'riskLastDayEnd','Risk final portfolio','rf','';...
         'placedInAnyPreferenceIniEnd','Placed to preference','ap','',1;...
    %     'placedInAddedIniEnd','Placed to added preference','apad','';...
    %     'placedInOldIniEnd','Placed to old preference','apold','';...
    %     'enrolledInAssigned','Enrolled in placed','ea','',1;...
    %     'enrolledInAssignedCOA','Enrolled in placed$|$placed','eaadcp','',1;...
    %     'declineOffer','Decline offer','do','',1;...
    %     'enrolledInAddedIniEnd','Enrolled in added','eaad','';...
    %     'addAnyUncongestedIniEnd','Add any uncongested','aaUncong','D. Congestion-related outcomes',1;...
    %     'placedInAnyPreferenceNWIE','Placed to uncong.','apUncong','',1;...
    %     'placedInAddedNoWaitlistIniEnd','Placed to uncong. added pref.','apadUncong','',1;...
    %     'enrolledInAddedNoWaitlistIniEnd','Enrolled in uncong. added','eaadUncong','',1;...
    %     'addAnyNoRiskIniEnd','Add any zero risk','aaNR','D2. Risk-related outcomes',1;...
    %     'placedInAnyPreferenceNoRiskIE','Placed to zero risk','apNR','',1;...
    %     'placedInAddedNoRiskIniEnd','Placed to zero risk added pref.','apadNR','',1;...
    %     'enrolledInAddedNoRiskIniEnd','Enrolled in zero risk added','eaadNR','',1;...
    %     'acceptOffer','Accept offer','ao','';...
    %     'preferenciaAsign_1','Order assigned','oa','';...
    };



if(anho==2020)
    outcomes=outcomes(not(ismember(outcomes(:,1),notIn2020)),:);
end



dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
for o=1:size(outcomes,1)
    
    dataRD.Properties.VariableDescriptions{outcomes{o,1}}=outcomes{o,2};
end


assert(allunique(outcomes(:,1)))
assert(allunique(outcomes(:,2)))
assert(allunique(outcomes(:,3)))

dataRD.riskPopupCentered=dataRD.riskPopup-dataRD.cutoff;
data=dataRD(:,[{'riskPopup','cutoff','riskPopupCentered'},outcomes(:,1)']);
inAnyPob=false(height(data),1);


% Generate stata commands
commands=cell(size(subpobs,1)*size(outcomes,1),2);
counter=1;


for p=1:size(subpobs,1)
    

    cutoffStr=subpobs{p,5};
    lbwStr=subpobs{p,6};
    rbwStr=subpobs{p,7};
     runVar=subpobs{p,8};
    
    data.(sprintf('subpop%i',p))=subpobs{p,2};
    inAnyPob=inAnyPob|subpobs{p,2};
    
    for o=1:size(outcomes,1)
        
        % Avoid computing outcomes that are not available:
        noCalcular=(strcmp(subpobs{p,3},'2020')|strcmp(subpobs{p,3},'2020comp'))&ismember(outcomes{o,1},notIn2020);
        
        if(not(noCalcular))
            
            
            % Full bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'full')||plotRDs)
                commands(counter,:)={sprintf('%s_%i_full',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, c(%s) p(2) h(%s %s)',outcomes{o,1},runVar,p,cutoffStr,lbwStr,rbwStr)};counter=counter+1;
            end
            % Full loccal fixed bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'local'))
                commands(counter,:)={sprintf('%s_%i_local',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, c(%s)  h(.1 .1)',outcomes{o,1},runVar,p,cutoffStr)};counter=counter+1;
            end
            % Full local auto bandwidht
            if(tableBandwidthComparison||strcmp(bandwidthTable,'localAuto'))
                commands(counter,:)={sprintf('%s_%i_localAuto',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, c(%s) ',outcomes{o,1},runVar,p,cutoffStr)};counter=counter+1;
            end
            
            % If with IV:
            if(outcomes{o,5}==1&&~isempty(subpobs{p,4}))
                
                endogenousVar=subpobs{p,4};
                
                % Full bandwidth
                if(strcmp(bandwidthTable,'full'))
                    commands(counter,:)={sprintf('%s_%i_full_iv',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, fuzzy(%s) c(%s) p(2) h(%s %s)',outcomes{o,1},runVar,p,endogenousVar,cutoffStr,lbwStr,rbwStr)};counter=counter+1;
                end
                % Full loccal fixed bandwidth
                if(strcmp(bandwidthTable,'local'))
                    commands(counter,:)={sprintf('%s_%i_local_iv',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, fuzzy(%s) c(%s)  h(.1 .1)',outcomes{o,1},runVar,p,endogenousVar,cutoffStr)};counter=counter+1;
                end
                % Full local auto bandwidht
                if(strcmp(bandwidthTable,'localAuto'))
                    commands(counter,:)={sprintf('%s_%i_localAuto_iv',outcomes{o,3},p),sprintf('rdrobust %s %s if subpop%i==1, fuzzy(%s) c(%s) ',outcomes{o,1},runVar,p,endogenousVar,cutoffStr)};counter=counter+1;
                end
            end
        end
    end
end
commands=commands(not(cellfun(@isempty,commands(:,1),'UniformOutput',true)),:);


if(false)
    % Run specific rd
    dataRD.auxRest=dataRD.restrAll&dataRD.mandatoryToAdd==0;
    a=stataCommand(sprintf('rdrobust lengthEnd riskPopup if auxRest==1, c(%s) p(1) h(.1 .1)',cutoffStr),dataRD)
    stataExpress(sprintf('rdrobust  DRiskExpostIniEnd riskPopup if auxRest==1, fuzzy(addAnyIniEnd) c(%s) p(1) h(.1 .1)',cutoffStr),data)
end

% Run Stata:
data=data(inAnyPob,:);
res=stataCommand(commands,data);


%% Make plots and table
close all

cantPobs=(size(subpobs,1));
cantRegs=(size(outcomes,1));

matrix_nl=nan(cantRegs,cantPobs);
matrix_nr=nan(cantRegs,cantPobs);
matrix_b=nan(cantRegs,cantPobs);
matrix_se=nan(cantRegs,cantPobs);

matrix_b_iv=nan(cantRegs,cantPobs);
matrix_se_iv=nan(cantRegs,cantPobs);
matrix_nl_iv=nan(cantRegs,cantPobs);
matrix_nr_iv=nan(cantRegs,cantPobs);

matrix_b_l=nan(cantRegs,cantPobs);
matrix_se_l=nan(cantRegs,cantPobs);

matrix_b_r=nan(cantRegs,cantPobs);
matrix_se_r=nan(cantRegs,cantPobs);

matrix_biasCorr_ci=nan(cantRegs,cantPobs,2);



regNames=fieldnames(res);

for p=1:cantPobs
    for r=1:cantRegs
        
        regName_i=sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable);
        if(ismember(regName_i,regNames))
            res_i=res.(regName_i);
            
            % Beta
            matrix_b(r,p)=res_i.tau_cl;
            % Sd Beta
            matrix_se(r,p)=res_i.se_tau_cl;
            
            % Beta left
            matrix_b_l(r,p)=res_i.beta_p_l(1);
            % Sd Beta left
            matrix_se_l(r,p)=sqrt(res_i.V_cl_l(1));
            
            % Beta right
            matrix_b_r(r,p)=res_i.beta_p_r(1);
            % Sd Beta right
            matrix_se_r(r,p)=sqrt(res_i.V_cl_r(1));
            
            % Bias corrected Confidence interval (alpha=5%)
            matrix_biasCorr_ci(r,p,1)=res_i.tau_bc-norminv(.975)*res_i.se_tau_rb;
            matrix_biasCorr_ci(r,p,2)=res_i.tau_bc+norminv(.975)*res_i.se_tau_rb;
            
            % N
            matrix_nl(r,p)=res_i.N_h_l;
            matrix_nr(r,p)=res_i.N_h_r;
            
            % If with IV:
            if(outcomes{r,5}==1&&~isempty(subpobs{p,4}))
                regName_i_iv=sprintf('%s_%i_%s_iv',outcomes{r,3},p,bandwidthTable);
                
                res_i_iv=res.(regName_i_iv);
                % Beta
                matrix_b_iv(r,p)=res_i_iv.tau_cl;
                % Sd Beta
                matrix_se_iv(r,p)=res_i_iv.se_tau_cl;
                
                % N
                matrix_nl_iv(r,p)=res_i_iv.N_h_l;
                matrix_nr_iv(r,p)=res_i_iv.N_h_r;
                
            end
            
            if(plotRDs)
                
                
                res_i=res.(sprintf('%s_%i_%s',outcomes{r,3},p,'full'));
                res_i_table=res.(sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable));
                
                plotsub=subpobs{p,2};
                
                figure
                plotRD(res_i,dataRD,plotsub,'otherPointEstimate',res_i_table);
                
                if(titlePlotRD)
                    title(sprintf('RD %s - %s',outcomes{r,2},subpobs{p,1}))
                    
                end
                
                if(savePlotRDs)
                    easyExport([dirPlots,sprintf(sprintf('%s_P_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                
                if(plotWithZoom)
                    % Plot with zoom
                    figure
                    plotRD(res_i,dataRD,plotsub,'newxlim',[.1,.5],'nb',100,'otherPointEstimate',res_i_table);
                    easyExport([dirPlots,sprintf(sprintf('%s_P_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
          
            else
                fprintf('Ojo: reg %s no existe para subpob %s\n',regName_i,subpobs{p,1});
            end
        end
    end
    
    %vector_Nl(1,p)=sum(dataRD.riskPopup<.3&subpobs{p,2});
    %vector_Nr(1,p)=sum(dataRD.riskPopup>.3&subpobs{p,2});
    
end





%% RD main table
matTable=nan(size(matrix_b).*[1 2]+[2 0]);
matTable_sd=nan(size(matrix_b).*[1 2]+[2 0]);

matTable(:,1:2:end)=[matrix_b;...
    matrix_nl(1,:)...
    ;matrix_nr(1,:)];

matTable(:,2:2:end)=[matrix_b_iv;...
    matrix_nl_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)...
    ;matrix_nr_iv(find(cell2mat(outcomes(:,5)),1,'first'),:)];


matTable_sd(:,1:2:end)=[matrix_se;nan(2,size(matrix_se,2))];
matTable_sd(:,2:2:end)=[matrix_se_iv;nan(2,size(matrix_se_iv,2))];

withInfo=not(all(isnan(matTable),1));


header=repmat({''},2,size(matTable,2));
header(1,1:2:end)=subpobs(:,1)';
header(1,2:2:end)=subpobs(:,1)';
header(2,2:2:end)={'IV'};


matTable=matTable(:,withInfo);
matTable_sd=matTable_sd(:,withInfo);
header=header(:,withInfo);

header=[{'Risky cutoff';''},header];

opt=struct;
opt.header=header;
opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=[outcomes(:,2);{'NL';'NR'}];
opt.addColumnNumber=true;

opt.filaFantasma=size(matTable,1)-2;
paneles=find(not(ismissing(outcomes(:,end-1))));
opt.panel=[num2cell(paneles-1),outcomes(paneles,end-1)];
%opt.alignmentFirstCol={'L{3cm}'};
opt.adjust=true;

opt.file=sprintf('%s_tablePopup%s',dirTable,anhoStr);

opt.title='RD Estimates of Platform Pop-Up Effects - 2017';
opt.label='tab:RDpopup2017';
opt.positionParameter='H';

opt.note=['Local linear RD estimates of pop-up effects from warning pop-up on application platform. Computed using triangular kernel with bandwidth 0.1. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. ',...
     'We report estimates in the pooled sample and for each different risky cutoff definition. IV (column 2) shows the instrumental variable specifications, where the endogenous regressor is the add any school indicator. ',...
     ''];

tabla=cell2latex(mat2cellstr(matTable,'precisionDecimal','%.3f','revisarFilas',true),'opts',opt);
if(compileLatexTable)
compileLatex(tabla)
end


